Binary black hole late inspiral: Simulations for gravitational wave observations 
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Coalescing binary black hole mergers are expected to be the strongest gravitational wave sources 
for ground-based interferometers, such as LIGO, VIRGO, and GEO600, as well as the space-based 
interferometer LISA. Until recently it has been impossible to reliably derive the predictions of 
General Relativity for the final merger stage, which takes place in the strong-field regime. Recent 
progress in numerical relativity simulations is, however, revolutionizing our understanding of these 
systems. We examine here the specific case of merging equal-mass Schwarzschild black holes in 
detail, presenting new simulations in which the black holes start in the late inspiral stage on orbits 
with very low eccentricity and evolve for ~ 1200Af through ~ 7 orbits before merging. We study 
the accuracy and consistency of our simulations and the resulting gravitational waveforms, which 
encompass ~ 14 cycles before merger, and highlight the importance of using frequency (rather 
than time) to set the physical reference when comparing models. Matching our results to PN 
calculations for the earlier parts of the inspiral provides a combined waveform with less than one 
cycle of accumulated phase error through the entire coalescence. Using this waveform, we calculate 
signal-to-noise ratios (SNRs) for iLIGO, adLIGO, and LISA, highlighting the contributions from the 
late-inspiral and merger-ringdown parts of the waveform, which can now be simulated numerically. 
Contour plots of SNR as a function of z and M show that adLIGO can achieve SNR > 10 for some 
intermediate-mass binary black holes (IMBBHs) out to z ~ 1, and that LISA can see massive binary 
black holes (MBBHs) in the range 3 x 10 4 < M/M Q < 10 7 at SNR > 100 out to the earliest epochs 
of structure formation at z > 15. 

PACS numbers: 04.25.Dm, 04.30.Db, 04.70.Bw, 04.80. Nn 95.30.Sf, 95.55.Ym 97.60.Lf 



I. INTRODUCTION 

Coalescing binary black holes (BBHs) are strong 
sources of gravitational radiation for ground-based de- 
tectors such as LIGO, VIRGO, and GEO600, and for the 
space-based LISA. This coalescence is generally consid- 
ered to proceed in three phases: inspiral, merger, and 
ringdown During the inspiral stage, the black holes 
are well separated and spiral together on quasicircular or- 
bits. This is followed by the merger stage, during which 
the black holes plunge towards the center and merge to- 
gether, forming a common horizon. This distorted rem- 
nant then 'rings down' to form a Kerr black hole by shed- 
ding its nonaxisymmetric modes as gravitational radia- 
tion. 

Different techniques are used to calculate the gravi- 
tational radiation from these stages. The inspiral can 
be treated analytically with post-Newtonian (PN) tech- 
niques; the gravitational waveform is a chirp, a sinusoid 
that increases in both amplitude and frequency. The 
later part of the ringdown can be handled analytically 
using black hole perturbation theory, with the result- 
ing waveforms being exponentially damped sinusoids of 
constant frequency. However, calculating the waveforms 



from the dynamical merger, which produces the highest 
luminosity signal, requires numerical relativity simula- 
tions of the full Einstein equations in three spatial di- 
mensions plus time. With the first generation of ground- 
based detectors now taking data and LISA moving for- 
ward, knowledge of the merger waveforms and their im- 
pact on detectability and parameter estimation is urgent. 

Recently there has been dramatic progress in numeri- 
cal relativity calculations of BBH mergers. The first full 
orbit of an equal mass nonspinning BBH was achieved 
nearly three years ago [2| (see also [3|), using a confor- 
mal formalism and comoving coordinates, with the black 
holes represented as punctures [H, and held fixed in the 
grid. This was followed about a year and a half later 
by the first simulation of the final orbit, merger, and 
ringdown [5]; this work was carried out using general- 
ized harmonic coordinates with excised black holes mov- 
ing through the computational domain 0, Q . Roughly 
six months later the moving puncture method, which is 
based on a conformal formalism and allows the puncture 
black holes to move across the grid without the need for 
excision [1,0], was introduced. This sim ple yet powerful 
technique proved to be highly effective [Tol [Til [T2I [l"3| , 
allowing the evolution of black holes through multiple 
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orbits, merger, and ringdown [3, EBj- Simulations with 
unequal masses [H, [It], EH and with spins [l!j[2(| quickly 
followed. Longer simulations with several orbits before 
merger have also been carried out using excision and har- 
monic techniques [2111, a llowing comparisons between the 
resulting waveforms [221 ]. 

This explosion of new work on BBH mergers has 
opened up exciting opportunities for applications of these 
results. A number of key applications center on gravita- 
tional wave data analysis. In particular, the need for 
accurate templates for the full gravitational wave sig- 
nal - encompassing the inspiral, merger, and ringdown 
phases - means it is essential to understand the relation- 
ship between PN waveforms and the results from numer- 
ical relativity [la, [2l|, [23| . And, with longer simulations 
now becoming available, it is becoming feasible to apply 
these waveforms to questions about the detectability of 
BBHs for various detectors [21]. In this paper, we focus 
especially on issues of accuracy in long simulations, start- 
ing in the late-inspiral regime, and their applications in 
gravitational wave data analysis. 

The most rapid advances cover the previously least un- 
derstood portion of the radiation, the final few cycles 
of radiation generated from near the "innermost stable 
circular orbit" (ISCO) and afterward, which we call the 
"merger-ringdown" . Already there has been considerable 
progress toward a full understanding of this important 
portion of the waveform through which the frequency 
sweeps by a factor of ~ 3 up to ringdown. A signifi- 
cant development in this regard was the demonstration of 
initial-data-independence of merger waveforms for equal- 
mass, spinless black holes [l5j . Today, there is general 
agreement that the merger of such a system produces a 
final Kerr BH remnant with spin a/m ~ 0.7, and that the 
amount of energy radiated in the form of gravitational 
waves, starting with the final few orbits and proceeding 
through the plunge, merger and ringdown, is -E ra d ~ 4% 
!, i, II EI El; see [24[ for a review. There is also 
consensus on the overall simple shape of the waveforms; 
detailed comparison of results among numerical relativ- 
ity groups has already begun [22j], with more to follow. 
The exploration of par ameter space, including various 
mass ratios [H, ITU Il8l| and spins [lj| |2(j, is now un- 
derway. Future efforts will involve pushing this frontier 
to increasingly complex mass-ratio and spin-orientation 
combinations, and to establishing initial-data indepen- 
dence across these parameters. 

Though more challenging, it is now becoming practical 
to simulate BBHs starting in the late inspiral as well. We 
report on new simulations covering roughly an additional 
factor of 3 in frequency before the "merger-ringdown" . 
We have simulated the coalescence of equal-mass, non- 
spinning black holes, starting in the late-inspiral regime, 
and continuing through the merger and ringdown. The 
black holes start on orbits with very low eccentricity and 
spiral together through ~ 7 orbits before merging. The 
resulting gravitational waveform has ~ 14 cycles, en- 
abling a robust consistency test with late inspiral PN 



waveforms as reported in (23|. In Sec. |II]we present an 
overview of our numerical results. Readers with an in- 
terest in the details of the simulations will find them in 
Sec. IIII1 where we discuss our methodology and the accu- 
racy of our simulations, as well as introduce the impor- 
tant notion of considering physical quantities as functions 
of frequency rather than time when comparing predic- 
tions from different computations. 

The rapid advances in understanding merger-ringdown 
waveforms with emerging simulation results for the late 
inspiral create a new context for considering gravitational 
wave observations of this part of the signal from such a 
system. In Sec. IIVI we revisit an analysis of merger de- 
tectability originally carried out in Flanagan and Hughes 
[l[ , now in the context of presentday and emerging wave- 
form modeling capabilities. In particular, we exam- 
ine the relative contributions to the signal-to-noise ratio 
(SNR) for the merger-ringdown and late-inspiral portions 
of equal-mass coalescence signals, as observed by initial 
LIGO (iLIGO), advanced LIGO (adLIGO), or LISA. We 
also plot contours of SNR as functions of redshift z and 
total mass for adLIGO and LISA, showing their abil- 
ity to detect BBHs in the cosmos. These SNR calcu- 
lations draw on a refined waveform model based on a 
best-estimate waveform using both numerical relativity 
and PN results. We conclude in Sec. [V] with a summary 
and discussion of our main results. 



II. OVERVIEW OF SIMULATION RESULTS 

Late-inspiral evolutions, covering more than a few or- 
bits prior to ringdown, are more challenging than shorter 
merger-ringdown simulations. In this regime there is a 
stronger requirement for numerical stability and a greater 
demand for computational resources. In addition, better 
accuracy is required to control the accumulated phase 
error, which in turn constrains the numerical error that 
can be tolerated in the rate of energy loss. 

Using the moving puncture technique [1, H E(3] , with 
the gauge evolution given by Eq.(17) and Eq.(26) in 
[l(]| |. we have simulated the evolution of a nonspinning 
equal-mass BBH starting at relatively wide separation, 
~ 1200A/ or ~ 7 orbits before the formation of a common 
event horizon. Here M is the total mass the system would 
have had when the black holes were very far apart and be- 
fore radiative losses became significant. We used fourth- 
order-accurate finite differencing techniques with adap- 
tive mesh refinement (AMR) to resolve the dynamics near 
the black holes and in the wave propagation region. We 
carried out three runs using similar grid refinement struc- 
tures, but at different resolutions: low (hf = 3M/64), 
medium (h f = 3M/80) and high (hf = M/32). Here, h f 
is the grid spacing in the regions with the highest reso- 
lution in each simulation, those being the regions around 
each black hole. Overall, the Hamiltonian constraint con- 
verges at fourth order, and the momentum constraint at 
better than second order, throughout the runs. 
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FIG. 1: Gravitational strain waveforms from the merger of equal-mass Schwarzschild black holes. The late part of the merger 
(t > — 50A/) is robustly determined and relatively easily calculable, while simulations of the late inspiral (early part of the 
waveform) are rapidly approaching the phasing accuracy required for observational applications [23|. The solid blue line shows 
our current "best" numerical waveform. The dashed red line shows a comparison waveform from a run starting with the same 
initial data as R4 in Ref. [l5| and the dash-dotted green curve shows the results from the highest resolution Rl run in Ref. (l5| |. 
All waveforms have been extracted at J? cxt = 40M and shifted in time so that the moment of maximum 7/14 radiation amplitude 
occurs at time t = 0. The initial coordinate distance between the punctures, di, is indicated in all cases. 



Fig. [T] shows the resulting gravitational waveform in 
the context of recent developments in black hole evolu- 
tions. The dashed (blue) curve gives the gravitational 
wave strain from the dominant I = 2, m = 2 spin- 
weighted spherical harmonic from our highest resolution 
simulation, extracted at i? cx t = 40M. This represents an 
observation made on the equatorial plane of the system, 
where only a single polarization component contributes 
to the measured strain. Time t = is taken to be the 
moment of peak radiation amplitude as measured by the 
Weyl curvature tp4] the symbols along the time axis mark 
the points at which the system reaches the ISCO calcu- 
lated for a Schwarzschild black hole (circle), EOB 3PN 
dHH (diamond), and adiabatic 3 PN [HElI (square). 
For comparison, we show two other waveforms from ear- 
lier simulations carried out by this group; both were ex- 
tracted at .R ex t = 40M and have been shifted in time 
and initial phase so (as in [151 ]) that the moment of peak 
ip/± radiation amplitude occurs at t = 0. As these dif- 
ferent simulations may radiate different fractions of the 
initial mass, we choose the mass to/ of the post-merger 
Kerr hole as the natural length scale for comparison (see 
discussion of Table U for details). Because of radiative 
losses, to/ 0.95M. 

The solid (red) curve shows the results from a model 
that starts ~ 550Af before merger with the same initial 
data as run R4 in Ref. [l5[ but using a different gauge. 
Note that the gauge conditions used for the dashed (blue) 
curve and the solid (red) curve are equivalent to those 
given by case #8 in Ref. [Toj , while the gauge conditions 
used in Ref. [HI] are equivalent to those given by case #3 
in Ref. [13] ■ The dot-dash (green) curve shows a simula- 
tion that starts ~ 200M before merger; this is the high 



resolution run Rl from Ref. |l5j. All three waveforms 
agree to within *~ 1% for the merger-ringdown burst part 
of the waveform, starting near t ~ — 50 Af. 

Astrophysically, BBHs in this relatively late stage of 
their evolution are expected to be traveling on nearly 
circular orbits, as any initial eccentricity would have 
been radiated away early in their evolution. In this new 
model, we start the black holes on nearly circular orbits 
with very small eccentricity e < 0.01, as defined below. 
The resulting black hole trajectories are shown in Fig. [21 
where the tracks mark the paths of the punctures; for 
clarity, only a single black hole from each simulation is 
shown. The dashed (blue) curve shows the results from 
our new high resolution long run; note that the early 
orbits are very nearly circular. The solid (red) curve 
shows the trajectory of the moderately long comparison 
run (shown by the dashed line in Fig. [I}. At early times, 
this model is significantly more eccentric than our new 
results. At later times this trajectory locks on to that of 
our new run, ~ 2.5 orbits before merger. 

The black hole separation as a function of time is 
shown for these two models in Fig. [31 The greater ec- 
centricity of the previous model (solid curve) is clearly 
distinguished here. We also show all three resolutions of 
our newest model. The slight deviations from the over- 
all smooth trend give an indication of the small amount 
of eccentricity in these latter simulations. Note however 
the differences between these three resolutions, particu- 
larly in the total time between the start of the simulation 
and the merger. Although significant, the differences in 
merger time appear to converge at a rate consistent with 
fourth-order error. 

In the next section we consider numerical techniques 
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FIG. 2: The trajectory of one of the binary system's black 
holes through ~ 7 revolutions before coalescence for our high 
resolution case is shown by the dotted line. The solid line 
gives the trajectory of the moderately long comparison run. 
The initial coordinate distance between the punctures, ck, is 
indicated in both cases. 
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FIG. 3: The coordinate separation between the puncture 
black holes is shown as a function of time. The solid line shows 
the results for the comparison run, which has relatively large 
eccentricity. The other curves show the three resolutions for 
our new simulations, all having noticeably less eccentricity. 
Note that equivalent gauge evolution equations were used in 
all four cases. 



and the fidelity of the simulations in more detail, includ- 
ing the convergence and conservation properties. We also 
carry out a detailed analysis of the resulting gravitational 
waveforms, underscoring the differences between using 
time and frequency to set references for comparing mod- 
els. Readers more interested in issues of detectability and 
SNR analyses for iLIGO, adLIGO, and LISA are invited 



III. SIMULATION DETAILS 



A. Numerical Methodology 

For initial data we used Brandt-Briigmann puncture 
data[3, generated by a second-order-accurate multigrid 
solver, AMRMG[28j. The puncture parameters were de- 
termined by the Tichy-Briigmann prescription for qua- 
sicircular initial data[29(, adjusted slightly, as informed 
by previous empirical experience, to reduce eccentricity. 
Our adjustment was simply to reduce the initial coordi- 
nate separation by roughly 2% while increasing the ini- 
tial momentum of each puncture such that the product 
of the initial momentum with the initial coordinate sep- 
aration remains constant. The success of this approach 
in giving a circular inspiral is roughly indicated by the 
puncture track (dotted line) in Fig. [2] and the curves in 
Fig. [3J The puncture track was computed by integrat- 
ing the equation x punc = —0(x punc ), where x punc is the 
position of the puncture, and (3 the shift. 

This data was evolved using standard Baumgarte- 
Shapiro-Shibata-Nakamura (BSSN) evolution equations, 
with the addition of dissipation terms as in [30( and 
constraint-damping terms as in [3l| in order to ensure 
robust stability. The dissipation terms, however, were 
tapered with Gaussian functions so as to vanish at each 
puncture position; this modification proved necessary 
for accuracy. The gauge condition was that recom- 
mended in [10l | for moving punctures. Time integration 
was performed with a fourth-order Runge-Kutta algo- 
rithm, and spatial derivatives with fourth-order-accurate 
finite-differencing stencils. For the outer boundary we 
employed a second-order-accurate Sommerfeld condition, 
pushed to \xi\ = 1536M to keep reflections far from the 
source. AMR was implemented via the software package 
PARAMESH [32|, HH , and interpolation between refine- 
ment regions was fifth-order-accurate. Note that we use 
AMR only to resolve the sources, and the mesh will pro- 
gressively become coarser far away from the sources. Al- 
though the radiation which reaches the outer boundary 
during the course of the simulation, with wavelengths of 
> 100M, will not be well resolved in this lowest refine- 
ment region of grid-spacing h — 32Af (in the highest 
resolution run), reflections from there are causally dis- 
connected from our extraction radii at R < 100AZ . We 
do not use AMR to follow the radiation with fine mesh; 
instead we require only that the fixed mesh resolution in 
the region of the extraction surfaces be sufficient to re- 
solve the waves there. For example the extraction surface 
at R — 60M, in our highest resolution simulation, spans 
regions with grid spacing h — IM and h — 2M. 
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FIG. 4: Convergence plot for the Hamiltonian constraint Ch- 
The top panel shows results from the finest grid and has been 
scaled so that for 2.5 order convergence the curves should 
superpose. The bottom panel shows results from the second 
finest grid and has been scaled so that for fourth order conver- 
gence the curves should superpose; the curves indeed appear 
to be fourth order convergent. 
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FIG. 5: Convergence plot for the Momentum constraint Cm- 
The top panel shows results from the finest grid and has been 
scaled so that for 2.5 order convergence the curves should su- 
perpose. The bottom panel shows results from the second 
finest grid and has been scaled so that for fourth-order con- 
vergence the curves should superpose; the curves appear less 
than fourth-order convergent but better than second-order 
convergent. 



B. Simulation Analysis 

The accuracy of the simulations was assessed by var- 
ious means. We first considered the L\ norm of the 
constraints in each refinement region, the grid structure 
having been designed to be commensurate for all resolu- 
tions; results for the finest (top panel) and second finest 
(bottom panel) refinement regions are plotted in Fig. [3] 
and Fig. [5] The Hamiltonian and momentum constraints 
were found to be convergent at an apparent order of 2.5 
in the finest grid, where the error at the puncture can be 
expected to dominate and fourth-order finite differenc- 
ing must break down due to the irregularity there. In all 
coarser regions the Hamiltonian constraint appears to be 
fourth-order-convergent, while the momentum constraint 
appears closer to second-order convergent. 

From each simulation we have measured the radia- 
tion in the form of the complex Weyl tensor component 
-04, specified consistently with [36| to leading order in 
1/r. The gravitational wave strain is related to 04 by 
— h+ + ih x = 2^4, and can be recovered by integra- 
tion, with the two complex integration constants cho- 
sen to keep the strain close to oscillating about zero. 
For some applications we also examine waveforms in 
the form of the strain rate v = h + + ih Xl the quan- 
tity from which radiation energy is directly obtained. 
To extract the waveform information from the simula- 
tion we define a series of coordinate spheres of differ- 
ent radii R ext /M G {40,60,100} on which we measure 
spin-weighted spherical harmonic components of ifii via 
a second-order algorithm described in [H, H3|. In this 
paper we focus exclusively on the I = 2, to = 2 compo- 
nent of the radiation, which mirrors the I = 2, to = — 2 
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FIG. 6: 04 waveform calculated at three different extraction 
radii, and scaled by the approximate Schwarzschild areal ra- 
dius. Times have been shifted according to the approximate 
Schwarzschild "tortoise coordinate" [34], |35J , as appropriate to 
compare the R ex t = 40M and R ex t = 60M waveforms with 
the R ext = 100M waveform. 



component because of equatorial symmetry. Other com- 
ponents are considerably smaller, containing ~ 1% as 
much energy (l5j . 

Fig. [5] compares waveforms from our high-resolution 
simulation extracted on each of our three extraction 
spheres with the i? cx t = 40M and i? G xt = 60M wave- 
forms shifted in time by the intervening propagation time 
derived based on a Schwarzschild black hole background. 
The generally good agreement for all three waveforms in- 
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TABLE I: The parameters of the final black hole (mass and 
spin parameter) formed by merger. Here m r and (a/m) r are 
calculated from the complex ringdown frequency cu r , while 
and (a/m)f are calculated from the radiation measured 
at the given radius. Analysis of these quantities at i? ex t = 
40A/ across the three different resolutions suggests that they 
are at least second-order convergent, and comparison across 
extraction radii gives error bars for the hf = M/32 results 
for m r , (a/m) r , rrif and (a/m)f of 1%, 1%, 0.3%, and 3%, 
respectively. 



hf 


i?cxt 


ringdown 
lu t m r (a/m) r 


conservation 
m f ( a / m )f 


3M/64 


40M 


0.551 - 0.083; 0.971 0.71 


0.953 0.70 


3M/80 


40M 
60M 
100M 


0.553 - 0.087; 0.945 0.68 
0.553 - 0.088; 0.931 0.66 
0.553 - 0.088i 0.931 0.66 


0.954 0.72 
0.955 0.70 
0.959 0.71 


M/32 


40M 
60M 
100M 


0.553 - 0.085; 0.953 0.69 
0.553 - 0.085; 0.953 0.69 
0.553 - 0.086; 0.945 0.68 


0.955 0.73 
0.955 0.71 
0.958 0.71 



dicates that potentially worrying subtleties related to the 
relatively close distance of the extraction spheres, such 
as nonlinear propagation effects, or tetrad-specification 
sensitivity, do not seriously affect the waveforms. On the 
other hand, some small differences are evident. For the 
early portion of the waveforms, shown in the inset, the 
results from the closest extraction sphere show slight dif- 
ferences in amplitude and phase, suggesting that 1 / i?^ xt 
radiation details are not yet insignificant here; this is not 
surprising given that the extraction radius in this case 
is only ~ 1/4 of a wavelength. On the other hand, the 
waveform from the farthest extraction radius shows signs 
of dissipation for the later higher-frequency portion of the 
waveform. This is because the radiation has propagated 
significantly farther, through a lower-resolution region on 
the computational grid, by the time it is measured at 
Rcxt = 100M. The resolution in most of the intermedi- 
ate region is h = 2M, about six points per cycle for the 
ringdown radiation. For the rest of the paper, we primar- 
ily employ the waveforms extracted at the intermediate 
distance i? ox t = 60M, which is only weakly affected by 
either of the above short- and long-wavelength effects. 

Different information in the waveforms provides inde- 
pendent ways to deduce the energy and angular momen- 
tum of the final black hole produced by the merger. The 
ringdown dynamics of the black hole after merger con- 
tains direct information about the mass and spin param- 
eter of the merged black hole, so that we can deduce m r 
and (a/m) r by measuring the frequency and decay rate 
of the ringdown radiation. Alternatively, we can mea- 
sure the radiation energy E ra d and angular momentum 
content J ra d at a given radius. Then, since we know the 
initial values, conservation of energy and angular mo- 
mentum imply the values of m/ = MadmU=o — E rac i 
and (a/m)f = {Jadm\i=o - Jrad)/m 2 f . Results from 



both methods are shown in Table fl] By comparing m/ 
with m r and (a/m) r with (a/m)f, we can verify conser- 
vation of energy and angular momentum in the simula- 
tions. In Table Q] we see that at the highest resolution, 
energy is conserved to within ~ 1% and angular momen- 
tum to within ~ 6%; and the best conservation is seen at 
Rext — 60M, where energy is conserved to within ~ 0.2% 
and angular momentum to within ~ 3%. 

To an excellent approximation, the I — 2, m — 2 har- 
monic of the radiation is polarized in the form expected 
for radiation generated by circular motion. The polar 
representation of the I = 2, m = 2 component of the 
complex waveform, ^4,22 = A^it) exp(i0^,(f)), is partic- 
ularly natural for circularly polarized radiation, for which 
A^p and to = dej)^ jdt vary only slowly. The angular po- 
larization frequency u then provides a meaningful instan- 
taneous frequency obtained directly from the radiation, 
corresponding to twice the angular frequency of orbital 
motion when the black holes are still separate. Because 
the radiation is measured in the weak-field region of our 
simulations, where gauge dependence is minimal, this 
polarization frequency provides gauge-invariant informa- 
tion about the binary dynamics. 

If the orbital motion is eccentric, this will leave an im- 
print in the radiation, causing a slight decrease in the 
polarization frequency of radiation generated near apo- 
center. We can recognize eccentric motion by identifying 
periodic deviations from a smooth monotonic trend in 
the time development of the polarization frequency ui{i). 

Specifically, we looked at the polarization frequency 
ui(t) = dcj)/dt calculated from the strain rate v(t) = 
V(t) exp(i(j)(t)) . We see generally similar results whether 
we use strain, strain rate, or ^4 to define the frequency, 
but specifically show the strain rate because it comes out 
smoother than -04 with respect to small waveform noise, 
but without noticeable detrending issues as in the strain. 
We fit the time dependence of the frequency curves to 
a quartic trend, u> c , plus an eccentric modulation of the 
form duj = Asin($o + ^o(* — *o) + &(t — ^o) 2 ) where 
the quantities A, <I>o, ^0 and (l are fitting constants 
and to — 61 AT is a time offset approximately account- 
ing for the propagation time to i? cx t = 60M. Ignoring 
the early part of the simulation where there are transient 
initial data effects, and the late part where the secular 
trend is very strong, we fit the curves over the time range 
250M < t < lOOOAf. We get similar results for eccen- 
tricity whether we apply high-frequency filtering to the 
waveform before fitting, though we show only low-pass 
filtered curves in Fig. [7l The results of the fitting are 
summarized in Table HI1 

Variations in the details of the fitting procedure, such 
as using strain instead of strain rate or smoothing high- 
frequency noise from u)(t), give results consistent to 
roughly the number of significant digits shown for A and 
Oo, though and $0 vary more significantly in some 
cases. The period of eccentric oscillations indicated by 
rio is about 1.5 times the initial orbital period, decreas- 
ing slightly at a rate comparable to the rate at which the 
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TABLE II: Values resulting from eccentricity fitting. The 
magnitude of the eccentricity in these simulations (as given 
by AM or eo) appears to be at least second-order convergent. 



Run 


AM 


n M 


ClM 2 


$o 


eo 


3M/64 


5 x 1(T 4 


0.017 


3 x 10" 6 


1 


0.005 


3M/80 


7 x 


0.016 


4 x 10~ B 


1 


0.007 


M/32 


8 x 10" 4 


0.017 


3 x 10" 6 


1 


0.008 




t/M 



FIG. 7: Angular frequencies with eccentricity removed, 
aligned such that the frequencies agree when the hf = M/32 
simulation is 1000M before the peak in 1^4. Also shown is the 
frequency from the hf — M/32 resolution before the eccen- 
tricity was removed. 



wave period grows. Allowing A to evolve in time did not 
result in clearly improved fits. 

From these fits we define eccentricity based on the ef- 
fect, in Keplerian dynamics, of small eccentricity on or- 
bital frequency, which should provide a good approxima- 
tion in the adiabatic regime. Kepler's second law (con- 
servation of angular momentum) implies that L oc r 2 uj 
is constant, from which it follows that the eccentric fre- 
quency deviation dto/uj is twice the eccentric radial de- 
viation dr/r, suggesting the unitless definition of eccen- 
tricity e = A/2oj c . Note that, to second order in e, this 
definition of eccentricity is also equivalent to that put 
forth in Ref. [38| in a post-Newtonian context. The con- 
stancy of A in our fit is consistent with a linear decrease 
in eccentricity as frequency increases. The initial eccen- 
tricity eo, obtained in this way, is also given in Table ITT1 

Even more interesting than the estimates for eccen- 
tricity provided this way are the residual curves u c (t) = 
Lu(t) — cko(t) of angular frequency vs. time with the ec- 
centric part subtracted out. Note that the strictly si- 
nusoidal nature of our fit to the eccentric modulations 
represented by cko(t) implies that the underlying secular 
trend should be preserved. Fig. [7] shows the results for 
each of our three simulations together with the unmod- 
ified frequency uj(t) for the high-resolution case. The 



plotted curves have also been filtered to remove some 
high-frequency simulation noise present in the early part 
of the simulations. The smooth, now monotonic, trend 
in the curves for u c (t) is an indication that most of the 
wiggles evident in uj(t) were consistent with our model 
for eccentric modulations duj{t), though the early part is 
affected by transient features related to the shortcomings 
of the initial data model. 

We can take uj c (t) as providing a record of the "harden- 
ing" process, as radiative losses bring the system through 
tightening orbits. The key effect of numerical simulation 
error evident in Fig. [7] is a slowing of the hardening rate 
at low resolutions, causing the final merger to be de- 
layed. This delay appears to converge at fourth order in 
resolution. Viewing the eccentric modulation as a small 
perturbation on the dynamics of an optimally noneccen- 
tric inspiral, we expect that these trends would also pro- 
vide a good approximation for the frequency evolution of 
a merger simulation begun with optimally noneccentric 
initial data. 



C. Waveform accuracy 

In this section we consider the accuracy of the wave- 
forms generated by our simulations, focusing primarily 
on the accuracy of waveform phasing information in the 
late-inspiral portion of our simulations. Over the course 
of many developmental simulations leading up to these 
results, we found that this early low-frequency part of 
the simulations is the most difficult to simulate accu- 
rately. This makes sense because the crucial dynamical 
details are in the slow loss of energy and angular mo- 
mentum to the relatively fine, evolving structure of the 
spacetime curvature, which ultimately comprises the ra- 
diation. Timescales are also longer for this part of the 
dynamics, requiring high accuracy over a large number 
of computational iterations. 

Fig. [8] compares the ■04 waveforms generated by our 
simulations at different resolutions. 

These waveforms have been aligned to coincide at the 
peak in ip^ for each of the simulations, as described in 
Ref. [IH . We show the waveform logarithmically because 
the amplitude changes by more than two orders of magni- 
tude through a run, part of what makes these simulations 
challenging. We are especially interested in the phase 
agreement among the different simulations, evident in the 
nearly vertical parts of the curves (which approach zero 
crossings). Aligned in this way, waveforms generated at 
all three resolutions are nearly exactly superposed after 
t ~ — 250 Af. The two higher-resolution simulations are 
nearly identical after t ~ — 500 Af and agree to within a 
small part of a cycle for the full portion of the simulation 
after the initial transient period in the first 100M of each 
run. On the other hand, the waveforms are easily distin- 
guished by large differences in the overall timing of each 
run, with the lowest resolution simulation going through 
a full additional orbit before merger. The high-frequency 
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FIG. 8: Waveform at three resolutions. They have been 
shifted in time and phase to agree at the peak of the wave. 
The agreement persists backwards in time, with the growing 
discrepancy in phase converging away with resolution. 




t/M 



FIG. 9: Strain rate phase. The high-resolution simulation 
goes through about 14 cycles before merger. The lower reso- 
lution simulations take longer to merge, and go through more 
cycles. 

noise evident in the first part of the simulations seems to 
be caused by reflections from our refinement boundary 
interfaces. 

We get a more direct view of phasing information by 
examining the waveforms in polar decomposition. In 
Fig. [9] we show the polarization phase derived from the 
strain rate waveform. 

This time we have aligned the simulations in time from 
the beginning, as should be appropriate for simulations 
of the same initial configuration. Later in the simula- 
tions, as the frequency increases the phase increases more 
rapidly, and the timing differences among the simulations 
lead to large phase differences. 

We quantitatively compare the phasing results in 
Fig. [TU1 showing the difference between the two higher- 
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FIG. 10: Strain rate phase three-point convergence. The 
higher resolution phase difference is shown with and without 
a phase shift to allow comparison of errors at similar dynam- 
ical points in the simulation. After shifting, the phase-error 
appears to be fourth-order convergent for much of the simu- 
lation. 



resolution runs compared with scaled versions of the dif- 
ference between the two lower resolution runs. We can 
now clearly see that the phase errors (measured this way) 
grow strongly in time. The lower-resolution difference 
has been rescaled two ways, such that they would be ex- 
pected to agree with the higher resolution difference for 
second- and fourth-order convergence, respectively. The 
comparison suggests phasing convergence somewhere be- 
tween second- and fourth-order over much of the run. 

Since the phase error grows strongly in time, it is in- 
teresting to consider the effect of timing errors in these 
convergence comparisons. The last curve in the plot 
shows the high-medium difference shifted in time by 
57M; that is, the high resolution and medium resolu- 
tion results have been differenced first, with their orig- 
inal time-dependence unaltered, and then the resulting 
difference has been shifted in time. 57M represents the 
approximate timing difference between the two higher 
resolution runs late in the simulation, as measured from 
the peaks in ip±. Viewed this way we see time- aligned 
phase differences taken from similar points in the phys- 
ical evolution of the medium resolution run, as repre- 
sented in the medium-low curve, and the high resolution 
run, as represented in the high-medium curve. As a re- 
sult, the runs appear to converge at a faster rate, closer 
to fourth-order. This unconventionally time-shifted plot 
is not intended to serve as a rigorous assessment of con- 
vergence, but is intended to illustrate that the timing 
differences can have a significant impact on convergence 
estimates. 

Above, we have compared our simulations made at dif- 
ferent resolutions by comparing the waveforms at equal 
points in time, with time aligned either at the beginning 
of the simulation (t = in the original run) or at the peak 
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FIG. 11: Frequency-based phase comparisons for runs at three 
resolutions. The solid curve represents the phase-difference 
between the hf = 3A//80 and hf = M/32 simulations; the 
dotted curve represents the phase-difference between the hf — 
3Af/64 and hf — 3A//80 simulations; the short-dashed curve 
represents the phase-difference between the hf — 3M/64 and 
hf — 3M/80 simulations, shifted vertically by a constant 
phase so as to agree with the higher resolution difference at 
LO c irif — 0.05423 (shown as vertical dot-dashed line); and the 
long-dashed curve shows a fit to u>~ 5 phasing error, which 
might be expected if there are energy conservation violations 
that are constant in time. The lower-resolution curves have 
been scaled so that they should superpose with the higher- 
resolution curves in the case of fourth-order convergence. The 
unshifted phase-difference curves appear better than fourth- 
order convergent, while shifting the lower-resolution curve 
makes the rate of convergence look closer to fourth-order. 



in ipi. Since errors cause the simulations to accelerate 
through their dynamical processes at somewhat different 
rates, such comparisons end up relating moments of quite 
different dynamics, and become less meaningful, depend- 
ing significantly on the reference time according to which 
the phases are compared. This sort of comparison would 
not be applicable at all when there is no clear way to 
physically align predictions in time, such as when com- 
paring numerical simulations with post-Newtonian (PN) 
results. We can avoid assigning a reference time by using 
the gravitational-wave frequency as the reference. 

For the case studied here, a quasicircular inspiral of 
comparable-mass black holes, it is appropriate to con- 
sider the waveform frequency u> — d(f>/dt as a gradually 
increasing monotonic function of time, ix> = u>(t). Though 
the actual waveform frequency of our simulations is not 
monotonic because of small eccentricities in the inspiral 
trajectories, we have shown that we can fit out the ec- 
centric deviations. The residual secular trend Lu c (t) pro- 
vides a monotonic frequency, which we can apply as an 
independent variable against which to compare various 
simulations. 

We show frequency-based phase differences among 
the simulations in Fig. 111! which allows us to com- 



pare the difference between the two higher-resolution 
simulations with that between the two lower-resolution 
simulations. If the simulation errors are fourth-order- 
convergent, then the low- medium difference should be ap- 
proximately 2.784 times the medium-high difference. As 
is evident from the medium-high "(3M/80-M/32)" and 
low-medium "(3M/64-3M/80)/2.784" curves in Fig. |TT| 
the errors appear slightly over-convergent with respect 
to fourth-order scaling. This may be caused by phase er- 
ror accrued early in the lowest-resolution (hf = 3M/64) 
simulation, due to difficulty in resolving high-frequency 
components in the spurious Bowen-York radiation, as 
well as in an initial gauge pulse, which dominate at this 
time. This lowest resolution may therefore not quite be 
in the convergent regime during this early part of the 
simulation. If the phases are adjusted by a constant such 
that they match at some point after the main part of the 
Bowen-York pulse has passed, then fourth-order scaling 
fits more closely. This is demonstrated in Fig. [11] by 
the "(3M/64-3M/80)/2.784 (shifted)" curve, which has 
been vertically phase-shifted by a constant so as to agree 
with the "(3M/80-M/32)" curve at urnif = 0.05423, the 
frequency lOOOAf before merger in our high- resolution 
simulation. 

As is also clear in the figure, phase error accumulates 
most significantly in the low- frequency portion of the sim- 
ulation, uifio < 0.08. This makes sense generally since 
the simulation spends much more time at lower frequen- 
cies. Consider, for instance, the effect of a small non- 
conservative leakage of energy SE from the simulation. 
When the black holes are well separated, the dynami- 
cal development manifested in the sweeping frequency is 
driven by a slow loss of energy and angular momentum. 
An energy leakage SE will change the frequency sweep- 
rate uj by Soj ~ (SE)/(dE/du>), where dE/dui indicates 
how the binding energy changes with frequency. Phase 
error 8<j) can be determined from the error in the sweep- 
rate by 

54* = 5 [-duj = - f^rSujduj (1) 
J w J u 

Applying the leading-PN-order expressions for dE/dcu 
and oj(uj) gives a result proportional to u>~ 5 SE. Indeed 
this dependence fits our phase differences rather well, as 
shown in Fig. QT] Note, however, that this does not sin- 
gle out energy nonconservation as the source of error, 
as other errors may produce similar effects. A similar 
leakage of angular momentum would lead to error pro- 
portional to oj~ 4 , which also fits reasonably well. 

The convergence evidenced by Fig. [11] suggests that 
we can apply Richardson extrapolation to estimate the 
difference of our high- resolution (A//32) run from the 
infinite-resolution limit. If we assumed the fourth-order 
convergence suggested by Fig. [Til the phase error esti- 
mate for this run would be 0.93 times the difference be- 
tween the phases from the 3A//80 and M/32 resolution 
runs. However we will simply take the more conservative 
estimate of the actual difference between these resolu- 
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tions. Note that during the last 1000M of our M/32 sim- 
ulation (i.e. from uimt — 0.05423 onwards), we estimate 
roughly two and a half radians of phase error accumulate, 
as measured with respect to frequency, which is less than 
half of a gravitational wave cycle. A benchmark for ac- 
cumulated waveform phasing errors is one-half of a cycle, 
because phase error exceeding this amount would lead to 
destructive interference in matched-filtering applications. 
For our high-resolution simulation, we estimate less than 
one-half cycle of gravitational wave phase error over the 
full simulated waveform, excluding the meaningless tran- 
sients in the first lOOAf. As in [23|, we estimate that 
these phasing errors are smaller than the implicit phas- 
ing difference between the 3PN and 3.5PN expansions of 
Cj{u)) after t ~ -300M (uj c M ~ 0.08). For our data anal- 
ysis considerations we will only be using the numerical 
waveform after this point, for which the estimated phase 
error is well below a half cycle. 

Frequency-based phase comparisons, such as we have 
presented here, are better suited than time-based phase 
differences, which depend strongly on where the wave- 
forms are chosen to be aligned in time. The relationship 
between the two can be understood by considering a one- 
parameter family of waveform results, with a parameter 
A representing model dependence, in this case the numer- 
ical grid-spacing. The waveforms would provide phase as 
a function of time 4>\(t), from which we can derive fre- 
quency u)\(t), which is monotonic for small eccentricity. 
Inverting to obtain t\(uj) one can derive the frequency- 
based phasing 4>\(uj) = (f>\(t\(u>)). Now considering vari- 
ations S = d/dX near A = 0, one finds the relationship 
between frequency- and time-based phase comparisons 

5(f>(ij) = 5<t>{t{uj)) + u5t{uo). (2) 

This sheds some light on the often confusing issue of 
time alignment in the time-based comparisons shown ear- 
lier. Specifically, for a waveform that sweeps significantly 
through frequency, time- and frequency-based phase dif- 
ferences will be most similar when the time-based phase 
differences are aligned so that 5t vanishes where to is 
largest. In our case, the net frequency-based phase dif- 
ferences in Fig. [TT] are closer to net phase differences with 
time aligned at the end of the waveform, as in Fig. [5J than 
when time is aligned at the lower-frequency beginning, as 
in Figs. 1- [IDI 

IV. APPLICATIONS TO GRAVITATIONAL 
WAVE OBSERVATIONS OF BBHS 

In the first part of this paper, we presented a state- 
of-the-art calculation of the late inspiral and merger of 
equal mass, nonspinning BBHs, starting ~ 7 orbits be- 
fore the peak radiation amplitude. In Ref. [23| we showed 
that there is a significant region over which the waveform 
from the M/32 numerical simulation presented above 
agrees with waveforms derived from PN calculations. In 
this section we will use the best of both treatments by 



joining the final segment of the BBH evolution with the 
PN approximation for the preceding inspiral, as shown 
in Fig. [TJ] We will show below that the phasing infor- 
mation for this waveform has an estimated accuracy of 
better than one-half of a gravitational wave cycle, making 
the waveform applicable in a variety of matched-filtering- 
based gravitational wave data analysis applications. We 
use this waveform to calculate the relative detectability 
of the inspiral and the merger-ringdown, as well as the 
detectability of the entire waveform, for iLIGO, adLIGO, 
and LISA. We show examples of characteristic signal 
strains for both classes of detectors. We also compute 
SNRs for astrophysically interesting BBH masses, high- 
lighting the importance of the late inspiral and merger 
parts of the signal. Modifications of our results that may 
arise from BBHs with spins and unequal masses are con- 
sidered. 



A. Waveform matching and SNR calculation 

The PN waveform used for all the analyses in this sec- 
tion is of order 3.5PN in phase (as derived in [3^, HoT ]) 
and 2.5PN beyond the quadrupole approximation in am- 
plitude (as derived in |41j). It should be noted that the 
PN approximation that we use for the phase is actually 
an expansion of the chirp rate u> in terms of the frequency 
ui, which we then numerically integrate. Direct numeri- 
cal integration of the 3.5PN expansion of the chirp rate, 
^3.5PJV, for example in the integrand cIu>(u)/uj3.5pn), does 
not strictly respect the PN approximation of the phase, 
as the latter would require additional 3.5PN expansion of 
the integrand itself. However, the phase obtained in this 
manner will have the same convergence properties as the 
original LU3.5PN expansion, which is arguably a more fun- 
damental PN quantity because of its close relationship to 
the rate of energy loss, E, from which it is derived. Ad- 
ditional expansion of the phase appears to compromise 
its accuracy. Using shorter runs it had been observed in 
[2~J | that the phase obtained from numerical integration 
of the PN expansion of the chirp rate seemed to agree 
better with numerical simulations than did the analytic 
expressions; more recently it was demonstrated in [23j . 
using runs extending well into the late inspiral and with 
the effects of eccentricity mitigated, that the numerically 
integrated phase appears to be converging to the numer- 
ical result at frequencies where the full 3.5PN expansion 
of the phase is clearly invalid. We therefore use the nu- 
merically integrated 3.5PN chirp rate for the phasing of 
the PN portion of our waveform. 

Fig.[12]shows our numerical waveform from SeclHlover- 
laid with the PN waveform that was just described. To 
generate a complete, mass-scalable waveform, we match 
the frequency of the numerical simulation to the PN pre- 
diction, adjusting the phases to also be equal at that 
point as shown in Fig. 1121 and connecting the two halves 
to make a single waveform. This is done by shifting the 
PN waveform until the frequency equals the numerically- 
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FIG. 12: Numerical simulation results for gravitational wave strain compared with post-Newtonian estimates. The waveform 
shown is from the high resolution numerical simulation presented in Sec. HT1 overlaid here with a PN waveform with 3.5PN order 
phasing and 2.5PN order amplitude accuracy. The combined waveform, joined at t = — 328A/ (circle) is applied in Sec. II VI to 
calculate signal-to-noise ratios for iLIGO, adLIGO, and LISA. 



predicted frequency at a time in the simulation where 
the accuracy of the numerical data first surpasses the 
accuracy of the PN approximation, as estimated in [23[. 
Specifically, [23[ predicts this point of equivalent accuracy 
to occur at Mw ~ 0.08, which corresponds to t = — 328M 
(shown by the circle in Fig. [T2|) . It is worth noting that 
there was no need to adjust the PN amplitude for con- 
tinuity. The amplitude agreement with the numerical 
simulation is so good, and hence the resulting amplitude 
is so nearly continuous, that the small discontinuity fails 
to produce any discernible artifacts in the Fourier trans- 
form h(f) of the resulting waveform. 

Having generated a waveform, it is informative to esti- 
mate the waveform's phasing accuracy over the course of 
the BBH evolution. Note in Fig. [TTJ that for the portion 
of our M/32 simulation that is used in the waveform, 
we estimate ~ 0.5 radians of phase error. If we take 
the difference between 3 and 3.5 PN terms to be an es- 
timate of the phase error as in j23j , we can assess the 
error for the PN portion of the waveform. It was shown 
in [12] that the analytic PN phase expression accumu- 
lates very little error, on the order of 0.1 radians, until 
Mtu ~ lx 10 . Beginning our numerical phase integra- 
tion at this point and evaluating up to Mcu — 0.08 yields 
a gravitational wave phase error of ~ 3.6 radians, such 
that the total accumulated phase error over the entire 
waveform is ~ 4 radians. As stated previously, an ac- 
cumulated waveform phasing error of less than 7r radians 
is the threshold below which wave-matching comparisons 
may be used for matched-filtering applications. We es- 
timate that our combined waveform meets this criterion 
after a frequency of about Mui ~ 0.01 up to the ringdown 
frequency, Mlo ~ 0.5. We therefore have a waveform with 
sufficient accuracy to be useful as a template for gravi- 
tational wave detection. While templates will ultimately 
be needed for cases of greater astrophysical interest, and 



still greater accuracy will be required for the template 
to be useful for the purpose of parameter estimation, the 
construction of this waveform illustrates that the field of 
numerical relativity has matured to the point of being 
capable of producing results that are useful for gravita- 
tional wave data analysis. 

The calculated waveform that we have just described 
is actually the total strain on the equatorial plane, where 
h x vanishes and therefore h + provides the only contri- 
bution. To get the optimally-oriented strain amplitude 
(which is the total strain passing an observer on the equa- 
torial axis), we multiply this result by 2y/2, which is the 
ratio of peak total gravitational wave amplitude to the 
amplitude of h + alone in the quadrupole approximation. 
We then simply divide by v5 in order to convert to an 
orientation-averaged waveform for our subsequent anal- 
yses. This factor can be understood by observing that 
orientation-averaging is fully equivalent to averaging over 
all sky positions of the detector from the perspective of 
the BBH, and such sky-averaging results in a factor 1 / v5 
in sensitivity [l[. 

The SNR is calculated assuming matched-filtering is 
performed on the data, and that the waveforms are per- 
fect copies of the embedded signal. In this case, the sky- 
and waveform-polarization-averaged SNR is given by 



< {SNR) 2 >= / d(ln/) 



^char \f) 



(3) 



is the characteristic signal 
(/) 



where h char (f) = 2/ h(f) 

strain and h n (f) = v5/i rms (/) = y/5fS n (f) is the rms 
of the detector noise fluctuations multiplied by V5 for 
sky-averaging, with h(f) and S n (f) being the Fourier 
transform of the signal strain and the power spectral den- 
sity of the detector noise, respectively [l[. 
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The waveform scales with luminosity distance Dl and 
total mass M as h^ar oc (1 + z)M/Dl, while the time 
axis for an observed wave, after redshifting, scales as 
t oc (1 + z)M, so that the waveform shown in Fig. [T2] 
is applicable over all total masses and redshifts. When 
needed we relate luminosity distance to redshift z using 
cosmological parameters consistent with the most recent 
Wilkinson Microwave Anisotropy Probe (WMAP) results 
(fi A = 1 - Q M = 0.72, h = 0.73) [H| and the relation 



D L (z) = 



(l + z)c 



dz' 



Ho Jo v/r! M (l + z') 3 + A 



(4) 



For cases where an impractically long time series would 
be needed to cover the band with an adequate sampling 
rate, the waveform is extended in Fourier space to still 
lower frequencies (and consequently back further in time) 
using the quadrupole formula, 



|ftquad(/)| 



1 



2vT5L> L 



[(1 + z)Mf 



1/6 



(5) 



The PN portion of the waveform continues slightly past 
where its Fourier transform deviates from the quadrupole 
expression for \h(f)\ by ~ 2%, at which point (O is used 
to extend \h(f)\ back as far as necessary. The PN seg- 
ment is truncated at a higher frequency than the lowest 
frequency component of its Fourier transform in order to 
eliminate edge effects. Finally, the quasinormal ringdown 
at the end of the numerical simulation is extended by fit- 
ting a damping coefficient and fundamental frequency to 
the data in order to mitigate edge effects at the high- 
frequency end of the Fourier transform. 



B. Observing Stellar BBHs and IMBBHs with 
iLIGO and adLIGO 

Ground-based interferometers are sensitive to rela- 
tively high frequency gravitational waves from coalesc- 
ing stellar mass (M < lO 2 Af0) and intermediate mass 
(IM) (1O 2 M < M < 10 3 Af Q ) BBHs. In this section, we 
apply our combined waveform to consider the response 
of iLIGO and adLIGO to BBH coalescence, illustrat- 
ing the importance of numerical simulation results for 
ground-based detectors. LIGO has facilities in Hanford, 
WA and Livingston, LA; each facility has an interfer- 
ometer consisting of two 4 km long arms, and Hanford 
also has a 2 km detector. Both the iLIGO and adLIGO 
detectors are designed to detect high-frequency gravi- 
tational waves, with iLIGO sensitive in the frequency 
range 40Hz < / < 8000Hz and adLIGO in the range 
14Hz < / < 10 3 Hz. Initial LIGO is currently operating 
at or near the initial design sensitivity in a year-long sci- 
entific data-taking run. Advanced LIGO is a planned 
upgrade that will increase the detector sensitivity by 
roughly an order of magnitude across the band. In addi- 
tion, adLIGO can be tuned to optimize its sensitivity for 
different sources. 
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FIG. 13: The iLIGO (dashed) and adLIGO (dash-dotted) 
rms noise amplitudes h n with the characteristic amplitudes 
h c har of 6 example sources (solid). The locations on each 
hchar corresponding to the peak ip4 amplitude (circle) and 1 
second before the peak in the observer's frame (filled circle), 
as well as t = — 50Af (square) and t = — 1000M (diamond) 
in the source's frame, are as marked. The mass given is the 
combined rest mass of each black hole. 



For our analysis of iLIGO, we used the design sensi- 
tivity to characterize the detector noise [44|. This sensi- 
tivity assumes that the noise is seismically limited below 
40 Hz, thermally limited between 40 and 150 Hz, and 
shot-limited above 150 Hz. For adLIGO, unlike iLIGO, 
we had a choice of tuning configurations. We used the 
wide-band tuning typically associated with burst sources 
because of its dramatically superior sensitivity at higher 
frequencies, where the merger portion from many sources 
is predicted to occur [45J . This yielded an improved SNR 
for most masses compared to tunings that were optimized 
for only the early inspiral portion of the coalescence. 

In Fig. [13l we show h c h ar for several sources plotted 
relative to the h rms sensitivity curves for iLIGO (dashed 
line) and adLIGO (dash-dotted line). We plot these val- 
ues because the height of h c h ar above h rms is an indicator 
of the SNR, as can be seen by inspecting Eq. [3l 

By rescaling we can calculate the SNR as a function 
of redshifted mass, and particular luminosity distance 
D L . In Fig. [H we plot the SNR achievable by iLIGO 
for sources at a luminosity distance Dl = 100 Mpc as a 
function of redshifted mass (1 + z)M. Here, the dashed 
line shows the SNR from the early inspiral in the time 
range — oo < t < — 1000M, which is roughly up to the 
start of our run. The dotted line shows the SNR for the 
late inspiral, -1000M < t < -50M, where t = -50M is 
approximately the time at which the merger burst begins. 
The thin solid line gives the SNR for -50M < t < oo, 
and encompasses the merger-ringdown part of the signal. 
The thick solid line shows the SNR from the entire wave- 
form. Note that the addition of the merger-ringdown 
waveforms increases the SNR and extends the detectable 
mass range significantly. The merger-ringdown portion 
t > — 50 Af dominates for all equal-mass nonspinning 
merger observations detectable with SNR larger than 10 
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(1+z) M 

FIG. 14: SNR for sources at luminosity distance Dl = 100 
Mpc plotted vs. redshifted mass for iLIGO. The contributions 
from -oo < t < -1000M (dashed), -1000 < t < -50 M 
(dotted), and —50M < t < oo (thinner solid), as well as the 
SNR from the entire waveform (thicker solid) are shown. 

at 100 Mpc. 

This type of plot was first made in Ref. [l[, and it 
is useful to compare our results with theirs. Our SNR 
calculations are based on a full waveform for the case of 
equal-mass, nonspinning black holes. The work in Ref. 
[l| was done before merger waveforms were calculated 
and thus is based on estimates for the merger-ringdown 
regime. For example, they estimated a merger radiation 
efficiency of ~ 10%, which is higher than our results but 
may well obtain for mergers with spin. Comparing their 
Fig. 4 for the SNR for iLIGO with our Fig. [14] we note 
that their curve for the inspiral includes the radiation up 
to the merger and so should be compared to the combi- 
nation (in quadrature) of our dashed and dotted curves. 
Our result for the merger SNR is somewhat smaller than 
theirs, due to the smaller amount of radiation emitted in 
our mergers. More recently, an analysis of SNR for iLIGO 
using numerical relativity waveforms for the merger and 
PN waveforms for the inspiral was made in Ref. [2l| ; our 
results in Fig. [H] are similar to what they report in their 
Fig. 22. 

Fig. [T5] shows the SNR for sources at Dl = 1 Gpc for 
adLIGO. Comparing with Fig. [Mj we see that adLIGO 
will have a significantly higher sensitivity to BBHs over 
iLIGO. This point is reinforced in Fig. Q1)J which shows 
contours of SNR for adLIGO as functions of redshift 
z and total mass M. We find that for M ~ 200M Q , 
adLIGO should be able to achieve an SNR greater than 
10 out to nearly z = 1 for equal-mass nonspinning bina- 
ries. From Fig. [TJ] it is evident that these high SNRs de- 
pend strongly on the merger-ringdown part of the wave- 
form t > -50M. 

It is important to note that astrophysical BBHs are 
likely to have mass ratios different from unity, and that 




(1+z) M 

FIG. 15: SNR for sources at luminosity distance Dl = 1 Gpc 
plotted vs. redshifted mass for adLIGO. The contributions 
from -oo < t < -1000M (dashed), -1000 < t < -50M 
(dotted), and — 50M < t < oo (thinner solid), as well as the 
SNR from the entire waveform (thicker solid) are shown. 
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FIG. 16: SNR contour plot with mass and redshift depen- 
dence for adLIGO. 



this will reduce the SNRs computed here for the equal- 
mass case. For stellar BBHs, current work 46j shows 
that the mass ratios are rather broadly distributed. The 
rates for such mergers may be low, ~ 2yr~ 1 for adLIGO, 
depending on the evolution of the original binary through 
the common envelope phase. For IMBBHs, mass ratios 
in the range 0.1 ^ m\jm% < 1 are expected to be the 
most relevant, with potential rates of ~ 10 per year [47l |. 
although these rates are far more uncertain than those 
for stellar BBHs. We can apply the mass scalings from 
Ref. [l| to show the effect of mass ratios on the computed 
SNRs; specifically, SNR - 77 1 / 2 for the inspiral, and SNR 
~ 1] for the merger and ringdown, where r\ = fi/M and 
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fi = m\vnilM is the reduced mass. Astrophysical BBHs 
are also expected to be spinning and this can potentially 
raise the SNR, for example if there is a spin-induced 
hangup that generates more gravitational wave cycles in 
the merger [19j. 



C. Observing MBBHs with LISA 

LISA, a proposed space-based interferometer consist- 
ing of three 5 million km long arms, will be sensi- 
tive to low-frequency gravitational waves from coalescing 
MBBHs in the band 3 x l(T 5 Hz < / < 1Hz. The coa- 
lescing massive binary black holes (MBBHs) that radiate 
in this band will have masses M > 10 4 M Q . 

Fig. [T7] shows hchar for several MBBHs plotted rela- 
tive to the LISA sensitivity curve. We used the "stan- 
dard" LISA sensitivity curve [H, for frequencies 
above 1 x 1CP 4 Hz, with shot and pointing noise con- 
tributions totaling 20pm/vTIz of laser phase noise. For 
3 x 10~ 5 Hz < / < 1 x 10~ 4 Hz, we employed a more 
conservative estimate of the acceleration noise than the 
one given in [48j], instead assuming a steeper amplitude 
spectral density that falls off as /~ 3 constrained to match 
the standard sensitivity curve at 1 x 10 _4 Hz [50]. Below 
3 x 10~ 5 Hz, we assume the detector has no sensitivity, 
which is a reflection of the uncertainty of the sensitivity 
at such low frequencies and our desire to make conserva- 
tive estimates. The sensitivity model assumes that there 
are no correlated noise sources, and its power character- 
ization corresponds to a round trip through one arm of 
the interferometer. 

MBBH sources can remain in-band for LISA over a 
very broad frequency range. Therefore, unlike the case 
of iLIGO and adLIGO, LISA sources nearly always re- 
quire the use of the quadrupole approximation procedure 
mentioned above to extend h c har to sufficiently low fre- 
quencies. Also, since more massive BBHs chirp more 
slowly, a MBBH could potentially be in LISA's sensitive 
band for much longer than the mission's lifetime. To 
prevent unrealistic SNR values due to this excessive in- 
tegration time, the quadrupole formula is only used to 
extend h c har to a low enough frequency such that the to- 
tal h c har used in our calculations corresponds to 3 years 
of data in the detector's frame, which is a conservative 
estimate of the expected mission lifetime. 

The SNR for LISA is shown as a function of redshiftcd 
mass, normalized for Dl = lOGpc in Fig. [TS] The bump 
in the curves is caused by the binary confusion noise. 
Again we see the enhancement of SNR from the merger- 
ringdown part (thin solid line) of the waveforms, and 
confirm the strikingly large values of SNR obtainable by 
LISA for these sources seen in [l[ and [2l[. For systems 
with redshiftcd mass (l + z)M <3xl0 4 , the early lnspi- 
ral t < — 1000M portion of the waveform dominates. The 
highest SNRs for equal-mass nonspinning mergers are ob- 
tained for systems with(l + z)M > 10 6 , again dominated 
by the merger-ringdown portion of the waveform. 
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FIG. 17: LISA rms noise amplitude h rms from the detector 
only (dashed) and from the detector combined with the antic- 
ipated white dwarf binary confusion (dash-dotted) [5lJ with 
the characteristic amplitudes h c har of three example sources 
(solid). The locations on each h c h ar curve corresponding to 
the peak ip4 amplitude (circle), 1 hour before the peak (filled 
circle), 1 day before the peak (circle with inscribed cross), 
and 1 month before the peak (circle with inscribed square) 
in the observer's frame, as well as t — — 5QM (square) and 
t — — 1000A/ (diamond) in the source's frame, are as marked. 
The mass given is the combined rest mass of each black hole. 
When necessary, the quadrupole approximation is used to ex- 
tend h c har backward in time 3 years before the peak ^4 am- 
plitude in the detector's frame (dotted). 

Contours of SNR for LISA are shown in Fig. [19] and 
demonstrate that LISA can observe MBBHs through- 
out the observable universe at large SNRs. We find 
it encouraging that, in addition to the large SNR val- 
ues predicted for LISA overall, some of the largest val- 
ues out to the largest redshifts occur in the mass range 
1O 5 M < M < 10 7 M Q where models of BBH popu- 
lations predict that the binaries can coalesce within a 
Hubble time [H3] and that the event rates for LISA are 
several per year [5^]. As discussed above, the effects of 
unequal masses will tend to decrease these SNR values, 
while spins may increase or decrease them. 

Even for nonoptimal configurations, the presence of 
an MBBH coalescence in the LISA data stream can 
dominate all the anticipated noise sources. Fig. l20l 
shows a simulation of LISA's response to the merger 
of equal-mass nonspinning black holes with total mass 
M = 10 5 M Q located at redshift z = 15, and oriented 
so that LISA lies in the system's equatorial plane, where 
the radiation is weakest. The SNR for a signal from such 
a source will be ~ 200, averaged over sky positions and 
polarizations (see Fig. [T9|). 

V. SUMMARY AND DISCUSSION 

Coalescing BBHs are expected to be the strongest 
sources for both ground-based interferometers as well 
as the space-based LISA. In particular, the strong-field 
merger portion of the gravitational wave signal produces 
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FIG. 18: SNR for sources at luminosity distance Dl = 10 
Gpc plotted vs. redshifted mass for LISA. The contributions 
from the early inspiral — oo < t < — 1000M (dashed), late 
inspiral —1000 < t < — 50A/ (dotted), and merger-ringdown 
— 50M < t < oo (thinner solid), as well as the SNR from the 
entire waveform (thicker solid) are shown. 
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FIG. 19: SNR contour plot with mass and redshift depen- 
dence for LISA. Note that MBBHs with masses M > 10 7 M Q 
may not coalesce within a Hubble time [52| . 

an intense burst of radiation and has the highest luminos- 
ity, emitting more energy per second than the combined 
starlight emitted in the observable universe. 

Recent breakthroughs in numerical relativity have 
opened a new era in understanding the late stages of bi- 
nary black hole coalescence. We now have a good under- 
standing of the merger-ringdown signal, starting ~ 50M 
before the peak radiation amplitude, for equal-mass non- 
spinning BBHs. The late inspiral evolution, that is, more 
than a few orbits prior to ringdown, is more challenging. 
Such simulations require better numerical stability and 
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FIG. 20: Simulated LISA data stream showing LISA's re- 
sponse to a system of two equal-mass black holes (M = 
10 5 Mq) located at redshift z=15 observed on the system's 
equatorial plane. The quantity plotted is an unequal arm 
Michelson interferometer observable "X" [54|. The LISA re- 
sponse and instrumental noise are realized using the LISA 
Simulator 55, 56], and colored noise was added to represent 
the unresolvable galactic binary foreground with the spectrum 
used in Ref . [511 ] . The inset shows the signal over a longer du- 
ration where low- frequency noise is evident. 



more computational resources, as well as higher accuracy 
to control the accumulated phase error. 

In this paper, we presented new simulations of equal- 
mass nonspinning BBHs starting in the late inspiral 
regime and covering approximately an additional fac- 
tor of three in frequency before the merger-ringdown. 
We carried out runs at three resolutions, hf — 
3M/64,3M/80, and M/32. Our runs start with rela- 
tively low eccentricity and show good convergence and 
conservation properties. We have demonstrated the sta- 
bility and accuracy of our simulations over the course of 
an unprecedented seven orbits. We also showed the value 
of using frequency (rather than time) to set a reference for 
the purpose of comparing results between runs as well as 
with the PN approximation. In recasting phase vs. fre- 
quency we have found particularly good agreement, not 
only between the runs but also with PN predictions. 

We have also matched the resulting gravitational wave- 
forms to PN calculations covering the earlier parts of the 
inspiral. The resulting full waveform has less than 3/4 
cycle of accumulated phase error over its entire frequency 
band. Using this waveform, we calculated the SNRs for 
iLIGO, adLIGO, and LISA. Our results confirm the im- 
portance of the merger-ringdown signal, which yields the 
highest values of SNR for the majority of equal-mass sig- 
nals |l|, |21j. We also show the SNR for the late inspi- 
ral regime, which numerical simulations are now begin- 
ning to address. The late inspiral dominates the SNR 
for iLIGO and adLIGO for the lower mass (< a few 
x 1OM ) stellar BBHs, and the SNR for LISA generated 
by MBBHs with M ~ 10 5 M Q . Contour plots of SNR as a 
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function of z and M show that adLIGO can achieve SNR 
> 10 for IMBBHs out to nearly z ~ 1, and that LISA 
can observe MBBHs at SNR > 100 out to the earliest 
epochs of structure formation at z > 15. 

Our work has focused on equal-mass nonspinning 
BBHs. Astrophysically, BBHs are expected to have un- 
equal masses and spins. In general, the effects of un- 
equal masses will tend to decrease the resulting SNRs, 
while spins can increase them. Calculations of merger- 
ringdown waveforms for several mass ratios, and for some 
spins, are currently available; we expect this to be a sig- 
nificant area of focus in the foreseeable future, both ex- 
panding the range of parameters studied and extending 
the duration of the resulting simulations. 

As computational technologies mature, simulations of 
the merger can be used in conjunction with gravitational 
wave observations to probe gravity in the arena of strong 
fields. In particular, if the binary masses and spins can 
be obtained from observations of the inspiral (as demon- 
strated, e.g., in Ref. [EtJ for LISA), numerical relativity 
can be used to calculate the merger waveform. This will 
allow a comparison between the predictions of general 
relativity - or indeed, any other theory of gravity used in 



a numerical simulation - with observations in the regime 
of very strong gravity. 
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